function [dof_map, V, T, u1, u2, err] = ex2(n_round, FE_Type, FE_Order, Young, nu, method)
% function [dof_map, V,T, u1, u2, err] = ex2(n_round, FE_Type, FE_Order, Young, nu, method)
%%
%  study the locking phonema where nu->0.5
%   the efficiency of FEM and spline method!
%
%
%  Test pass on 11, Aug, 2012
%
%  Author: Dr. Xian-Liang Hu
%

if nargin < 1
    FE_Type = 'BB';  % default FE basis.
    FE_Order = 3;   % default polynomial order.
end

if nargin < 3
    % default physical parameters for elastic problems
    Young = 250; 
    nu = 0.25;
    % lambda = 100;
    % mu = 100;
end
R = 1/2;

if nargin < 5
    method = 2;  % using FEM style by default
end

% %% for output dofs v.s. errors
% dof_count = zeros(10,1);
% err_L2 = zeros(10,1);
% err_inf = zeros(10,1);

Quad_Order = 13;

% generate initial mesh
load square_circle;
[T, E, ET, TE] = build_fem_mesh(V,T);
for round = 1:n_round
    [V, T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET);
    % we have to pull the boundary nodes on inner circle in this example
    bnd = find(ET(:,2) == 0);
    bnd_V = unique([E(bnd,1); E(bnd,2)]);
    x = V(bnd_V,1); y = V(bnd_V, 2);
    dist = sqrt(x.*x + y.*y);
    inner =  dist < 0.6;
    x = V(bnd_V(inner),1);
    k = V(bnd_V(inner),2)./x;
    x_new = (x*R)./(abs(x).*sqrt(1+k.*k));
    V(bnd_V(inner),1) = x_new;   % map the nodes onto the circle (0,0,R); 
    V(bnd_V(inner),2) = k.*x_new;
end
% ZERO=zeros(size(V,1),1); trisurf(T,V(:,1),V(:,2),ZERO);view(0,90);
fprintf('There are %d triangles.\n', size(T,1));

% n_elem = size(T,1); n_dof_per_elem = (FE_Order + 1)*(FE_Order + 2)/2;

%%%%%%%%
%  
bdr_Dirichlet = find(ET(:,2)==0);
bdr_Neumann = [];

[dof_map, u1, u2, t] = elastic_fem(V, T, TE, ET, FE_Type, FE_Order, Quad_Order, method, ...
                                     Young, nu, @fun_f, bdr_Dirichlet, bdr_Neumann, @fun_u, R);


%% ����������

%%% this is for the linear FE case:
% n_dof = max(max(dof_map)); 
% u = zeros(n_dof,1); u(dof_map) = u1;
% trisurf(T,V(:,1),V(:,2),u);

%%% order d finite element solutions:
disp_d = 3; tri_temp = template_mesh_tri(disp_d);
[u_h, Tris1, Points1] = fe_solution_bb(V,T, u1, FE_Order, tri_temp, disp_d);
[v_h, Tris2, Points2] = fe_solution_bb(V,T, u2, FE_Order, tri_temp, disp_d);

% analytic solutions
[uu, vv] = fun_u(Points1(:,1),Points1(:,2), Young, nu, R);

% cal the L_inf error
erru = max(max(abs(u_h - uu)));
errv = max(max(abs(v_h - vv)));

% return the maximum
err = max(erru, errv);

% % visualization and print information
% subplot(1,2,1); 
trisurf(Tris1, Points1(:,1), Points1(:,2), u_h);  %plot error for u
% subplot(1,2,2); trisurf(Tris1, Points1(:,1), Points1(:,2), v_h - vv);  %plot error for v


fprintf('The inf norm or error(u) = %e,  error(v) = %e.\n',  erru, errv);

end


%%%%%%%%%%%%%%%%%%%%%%%%
%  analytic solution, also used as boundary function
% 
function [u_1, u_2] = fun_u(xx, yy, E, nu, R)
   r2 = (xx.*xx + yy.*yy);
   u_1 = r2/100 + (1 - 1/100)*R*R;
   u_2 = u_1;
end


%%%%%%%%%%%%%%%%%%%%%%%%
%  The right hand side function
% 
function [f_1, f_2] = fun_f(xx, yy, E, nu, R)
    [n_row, n_col] = size(xx);
    f_1 = -8*ones(n_row, n_col);
    f_2 = f_1;
end